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A leading approach to improving the accuracy on numerical relativity simulations of black hole 
systems is through fixed or adaptive mesh refinement techniques. We describe a generic numeri- 
cal error which manifests as slowly converging, artificial reflections from refinement boundaries in 
a broad class of mesh-refinement implementations, potentially limiting the effectiveness of mesh- 
refinement techniques for some numerical relativity applications. We elucidate this numerical effect 
by presenting a model problem which exhibits the phenomenon, but which is simple enough that 
its numerical error can be understood analytically. Our analysis shows that the effect is caused by 
variations in finite differencing error generated across low and high resolution regions, and that its 
slow convergence is caused by the presence of dramatic speed differences among propagation modes 
typical of 3+1 relativity. Lastly, we resolve the problem, presenting a class of finite-differencing 
stencil modifications which eliminate this pathology in both our model problem and in numerical 
relativity examples. 

PACS numbers: put pacs numbers here 


I. INTRODUCTION 

Recent years have seen a dramatic rise in opportu- 
nities for observing strong-field gravitational dynamics. 
New observations of dense black-hole-like objects, both 
at stellar, and supermassive scales are increasingly fre- 
quent. Anticipated gravitational wave observations by 
ground-based and space-based detectors are expected to 
capture information about these objects at moments of 
the strongest gravitational interactions. Interpretation 
of data from any such observations will depend on theo- 
retical modeling of the strong-field interactions of dense 
black-hole-like objects in the process of generating grav- 
itational radiation. General Relativity is the standard 
model for describing gravitational interactions and wave 
generation. Unfortunately, the predictions of General 
Relativity for such cases are not yet fully understood, 
but will depend on 3-D numerical relativity computer 
simulations. 

While numerical relativity has progressed markedly in 
recent years, significant improvements in the fidelity of 
models for events such as binary black hole coalescence 
will be essential for the full interpretation of upcoming 
observations. There are many facets to the problem of 
improving such simulation, including optimally formu- 
lating of Einstein’s equations, properly handling bound- 
aries, handling black hole singularities, handling con- 
straint violations, and making judicious gauge choices. 
There are also basic numerical issues concerning how 
to, with finite resources, perform such high-fidelity 3-D 
simualtions with strong short wavelength gravitational 
features near the sources, and weak but critical long- 
wavelength gravitational (radiation) features emerging in 
a large domain. Approaches to this latter class of issues 
include, developing higher-order accurate finite differenc- 
ing methods, spectral methods, mesh-refinement tech- 
niques and other forms of numerical patching techniques. 

We focus here on resolving a limitation which has 
arisen in our work on mesh refinement techniques, but 


which may have analogues in other numerical patch- 
ing treatments as well. Mesh refinement techniques di- 
vide the computational domain into regions with sepa- 
rate computational grids which can be of higher resolu- 
tion in some regions than others. Such approaches in- 
vlove mesh-structure interfaces across which the details 
of the finite-differencing treatment suddently change. In- 
evitably, these interfaces contribute to computational er- 
ror, manifesting such effect as “reflections” off the in- 
terfaces. Considerable attention is given to implement- 
ing “clean” interfaces, which generate small error, com- 
pared to error generated in the bulk regions. At mini- 
mum, this requires interface induced error to converge at 
least as rapidly as the bulk error. For some classes of 
black hole evolutions, with non- vanishing shift advection 
terms, we have typically observed large error propagat- 
ing from mesh interfaces. This problem forms the focus 
of this paper. 


II. INTERFACE PERFORMANCE 

The pertinent features of our numerical scheme are 
as follows. We are solving the 3+1 BSSN formulation 
of Einstein’s equations. Our gauge condition is numer- 
ically determined, typically with some variation of the 
1+log slicing and hyperbolic Gamma-driver shift evolu- 
tion equations. All spatial derivatives are computed by 
2nd order accurate, centered differencing, save for ad- 
vection derivatives, for which we use 2nd order accurate 
upwinded differencing. 

To resolve any sources, such as black holes, adequately, 
whilst pushing the computational grid boundary suffi- 
ciently far away, we use fixed mesh refinement (FMR), 
as implemented by a software package for this purpose 
called Paramesh. With this implementation, the resolu- 
tions of two adjacent refinement regions always differ by 
a factor of two; i.e. the grid-spacing of the coarser region 
is double what it is in the finer region. ’’Ghostzones” or 
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FIG. 1: Convergence plot for the Hamiltonian constraint H at 
time t = 4 M. There is a single puncture black hole centered 
at the origin, and refinement boundaries at |x<| = 1 M, 2 M, 
and 4M. The finest grid spacing hf for each simulation is 
indicated in the figure, and the higher resolution is multiplied 
by a factor of 4, as appropriate for demonstrating 2nd order 
convergence. 
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where a = a — 1, p\ = p l - Pq, (with Pq assumed 
spatially uniform for simplicity), and hij = 7 ij - <5^. 

For a one-dimensional problem, if we assume plane- 
wave solutions (generalizable by Fourier analysis), then 
the above system of equations can be written in the form: 


’’guardcells”, typically two layers, are required to provide 
buffering between refinement levels. These guardcells are 
filled in by interpolation. 

Whether the simulation performs adequately in the 
presence of refinement interfaces is a question of partic- 
ular concern to us. Of course, refinement boundaries 
are sources for reflection, but these reflections are, gen- 
erally speaking, satisfactorily convergent and often neg- 
ligible. An exception has plagued us in the case of a 
non-negligible shift across a refinement boundary. In this 
case, we observe reflection that propagates at a velocity 
of - P % , that is not clearly converging at an appropriate 
rate, and that, in some situations, dominates the error. 
This pathology is most evident in the Hamiltonian con- 
straint as computed from the evolution of a single Brill- 
Lindquist puncture black hole (Fig. 1). 


III. A MODEL PROBLEM 


A. Linearized BSSN system 


To understand the source of the p - speed error exhib- 
ited in the last section, we begin by considering a lin- 
earized BSSN system with 1+log slicing and hyperbolic 
Gamma-driver shift. The system of equations is: 
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The eigenvalues are 0,~Po ,-fio + 1 ,~Po ~ ^,~\Po + 
WM +8, -^/?o - \\fp% +8,-§/? 0 + \y/Pl + 4, and 

— \Po ~ 3 VM + 4. 

If |A) is the eigenvector associated with eigenvalue A, 
and (A| is defined such that | A) (A| = 1, then 
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By substituting Eq. (12) into Eq. (9) the system of 
evolution equations can be decomposed into a series of 
advection terms, each associated with a characteristic ve- 
locity equal to one of the eigenvalues. Note that, assum- 
ing Po 1, the /3-speed mode is much slower than any of 
the other non- zero-speed modes. As we are particularly 
concerned with this mode, it is instructive to write the 
evolution equation thusly: 
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FIG. 2: Convergence plot for the error in ip. There is a re- 
finement boundary at x = 50. The coarse grid spacing h of 
each simulation is indicated in the figure, and the medium 
and high resolution curves have been multiplied by 4 and 16 
respectively, as appropriate to demonstrate 2nd order conver- 
gence. 


This equation indicates that disturbances originating 
in </> and T can propagate in (f> and h at /3-speed and fur- 
ther, that this is the only means of /3-speed propagation 
allowed by this system. Thus </> is uniquely significant in 
both generating and propagating /3-speed modes. 


B. ’’Bumpy” system 

Less formally it is clear that Eq. (13) resembles the <fi 
evolution equation in the BSSN system. The simplest 
model for the generation and propagation of the /3-speed 
modes would be one dimensional advection equation with 
an additional driving term, 

= P?^( x ,t) + (14) 

In our numerical simulations, it appears that the reflected 
’bumps’ may be triggered by a rapidly propagating gauge 
pulse which propagates outward early in the simulations. 
For our model problem, which we will term ’’Bumpy” 
because of its most salient feature, we will drive the ad- 
vection equation with a pulse propagating at speed u, sig- 
nificantly faster than the advection speed /3. We thus let 
/(x, t) = f(x — vt) where v is larger than /3 and both are 
assumed positive. This model equation is simple enough 
to be solved directly, 

t) = f f(x + (t - t')0 , t')dt f 
Jo 

= -ZL (F(x-vt)-F(x + pt)). (15) 

For definiteness we can take the driving term to be a 
Gaussian pulse, /(x) = exp(— x 2 ), implying F(x) = 
( V y/ ( 7r )/2)erf(x). For a localized pulse, such as this, the 


F(x + fit) term in the exact solution will be negligible in 
the region of interest, near x = xq. 

To test if this model problem is sufficient, we have nu- 
merically evolved Eq. (14) using a 1-D finite difference 
code with a resolutions dx = h for x < Xo and dx = 2 h 
for x > xo, realizing a refinement boundary at x = xo. 
We evolved over the domain 0 < x < 100 with peri- 
odic boundary conditions, using three-point upwind finite 
difference stencil and a mesh refinement scheme similar 
to that in our numerical simulations with PARAMESH. 
Blah Blah figure Blah Blah. ... not strongly affected by 
variations in the time-integration method or by changing 
amoung mesh-refinement interfacing schemes which are 
consistent with the overall second-order finite differenc- 
ing accuracy. 


IV. A FINITE DIFFERENCE ANALYSIS 

In this section we attempt to understand the numer- 
ical behavior of the Bumpy problem Eq. (14) by con- 
structing here an analytic model for the numerical error 
in our simulations. Noting that the time discretization, 
and the details of the interpolation scheme used in apply- 
ing refinement conditions seem not to be directly linked 
to the problematic error features we see, we model the 
error with as few assumptions as possible about these de- 
tails. We treat the numerical errors continuously in time, 
and we consider the effects of spatial finite differencing 
in terms of a continuous field (p(x,t,h) representing the 
numerical solution at (fine-grid) resolution dx = h. We 
then expand tp in orders of h, 

y>(x, t, h) = ip e (x , t) + h 2 (p 2 (oo, t) -fi 0(/i 3 ), (16) 

Where < p c (x,£) is understood to be the exact solution 
Eq. (15). We consider the effect of our finite differencing 
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scheme by replacing the spatial deriative appearing Jn 
Eq. (14) with a suitable finite difference operator Dk- 
Thus, 


x,t,h ) = /3D h <p(x,t,h) + f{x,t). (17) 

As above, we include a refinement jump in the grid at 
x = Xoj represented here by applying a coarser version 
of the finite difference stencil in the x > xo part of the 
spatial domain. This will be consistent with any refine- 
ment interface algorithm which applies the same finite 
difference stencil in both the coarse and fine regions, and 
applies a interpolative guard-cell filling alrogithm at the 
interfaces which leads to finite differences at the inter- 
face which are consistently second-order accurate. On a 
uniform grid the second-order error term of a finite first 
derivative is generally proportional to the third derivative 
of the field, 

Dk<p(x,t, h) = d x <p(x,t,h) 4- e2h 2 dl(p(x,t,h) 4- 0(h 3 ). 

For the specific upwind differencing operator used in Sec- 
tion III, the stencil of which is, 

Dk = iH +2£ ' , “5 £ “) (18) 

where E is the spatial translation operator defined such 
that Ehf{x) = f(x 4- h), the constant error coefficient 
turns out to be e<i — —1/3. Including the refinement 
jump, we have 

Dh<P = D h <p + &(x - x 0 )(D 2 h -D h )tp (19) 

— d x tp + (1 + 30(x - xo))e 2 h 2 dl<p + 0(h 3 ). 
Substituting into Eq. (17), and rearranging, yields 

h 2 (p 2 {x,t,h) = [-(p e (x,t) + 0d x <p e {x,t) + f(x,t)] 
+h 2 0[d x ip2(x, t, h) + 

e2(l + 30(z - xq )) dl<p e (x,t)] 
+0(h 3 ). (20) 

Then noting that the first term vanishes, and taking the 
limit h — ► 0, we derive 

02 OM) = 0d z <p 2 {x,t)+ (21) 

pe 2 (l 4* 30(x - x 0 ))£%tp e (x, t), 

where we have used the notation <£2 (x,t) = (p2(x,t,0). 
This is our model for generation and propagation of fi- 
nite differencing error in the numerical model problem in 
Section III. 

Since Eq. (21) is of the same form as Eq. (21) we can 
solve it the like fashion. We leave out the negligible F(x 4- 
0t) in <p e , substituting 

f{x,t) = /?e 2 (l 4- 30(x - xo))d%<p e (x,t) 

~ -J?L.(l + ZQ(x-xo))%F(x-vt) 
v 4- p 



into the first line of Eq. (15). Then, performing the in- 
tegral with careful attention to the presence of the step 
function, we get the solution 


<p 2 (x,t) = (s(x - vt) - s(x 4- (3t))(l 4- 3©(x - xq)) 


-i-3s(-~(x - x 0 4- P(t - ^)) x 
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xo y 
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As above, the s{x 4- /?£), thus the first term in Eq. (22) is 
effectively the differencing error associated with the up- 
sweep in (p e which propagates across the grid at speed v. 
This part grows four times larger in the coarse x > xo re- 
gion. The second term is quite interesting, it propagates 
in the reverse direction at speed /?. It has the same s(x) 
shape as the forward propagating component, but it is 
reversed, and contracted by a factor /?/zz. It is timed to 
orginate at the interface as the first pulse crosses. The 
two terms combine make the full solution continuous at 
the interface. Heuristically, one could say that the second 
term is caused by the discontinuity in the differencing er- 
ror, generated in order to produce a regular solution, and 
then, necessarily advecting away as required by the orig- 
inal model equation Eq. (14). It is contracted because it 
propagates at a different speed than the first term, but 
their time dependences must match at the interface. 


V. IMPROVING THE DIFFERENCING 
STENCILS 

The results of the last section clearly show that the re- 
flected error in the Bumpy system is overwhelmingly due 
to the discontinuity in the differencing operator. More 
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specifically, the reflection is related to the discontinuity 
in the 2nd order truncation error. This observation sug- 
gests a solution. 

Consider using modified differencing stencils such that 
the coefficient of the 2nd order truncation error in the fine 
grid region is multiplied by a constant factor qo while the 
coefficient of the 2nd order truncation error in the coarse 
grid region is multiplied by a constant factor qi. Then 
the previously constant coefficient e 2 of the last section 
becomes: 

e 2 -+ [go + (<?i - qo)0(x - X 0 )]e 2 (23) 

Repeating the steps of the last section for this new 
truncation error, the solution for the 2nd order error in 
< p becomes, 



FIG. 4: Comparison of the error in (p reflected from the re- 
finement boundary at x = 50 in the case of a 2nd order con- 
ventional stencil and a 2nd order MAD stencil. 


<P 2 (x,t) = (s(x - vt) - s(x + /3t))(q 0 + (4gi - qo)Q(x - x 0 )) 
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stencil combined with a 3rd order lopsided stencil, for ex- 
ample, the resulting stencil has the form: 


(©(x — xq + fit) ~ 0(x - x 0 )) 

Thus for the spatially blueshifted, reflected error, we now 
obtain, 


ifiref = (4gi - q 0 )s(-^(x - Xo + 0(t - ^))) (25) 

Clearly, this reflected error can be eliminated simply 
by choosing qo = qi = 0. This choice corresponds to 
using higher order stencils, say 3rd or 4th order accu- 
rate, throughout the entire grid. However, it has been 
observed that higher order differencing can some times 
create instabilities in a strong field region (which would 
typically correspond to the finest grid region). So it may 
be useful to consider an alternative solution. 

Note that with the choice qo ~ 1 and q\ = the 
reflected error vanishes. More generally, of course, any 
choice of qo and qi that makes the truncation error con- 
tinuous across the refinement boundary will remove the 
2nd order reflection. In particular, the choice 



will work in the nth refinement region of a grid with an 
arbitrary number of refinement levels, where ho is as- 
sumed to be the grid-spacing in the finest region. We 
call such a scheme mesh-adapted differencing (MAD). 

A 2nd order MAD stencil can be obtained simply by 
linearly combining a 2nd order accurate stencil with a 
higher order accurate stencil, as follows: 

D hn =q n DM+(l-q n )DW (27) 

where the superscripted numbers in square brackets rep- 
resent the order of accuracy of the differencing operator, 
and j > 2. In the particular case of a 2nd order upwinded 
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Implementing this MAD stencil in the Bumpy code, it 
is proven phenomenologically (and phrenologically) that 
the dominant reflected error is eliminated (Fig. 4). 


VI. RESULTS FOR BLACK HOLE 
EVOLUTIONS 

Returning to our Einstein solver, we have found that 
3rd or higher order accurate differencing of advection 
terms is unstable in the vicinity of black hole punctures, 
if the stencil has any points on the ’’downwind” side (op- 
posite the shift, towards the puncture). If the 3rd or 
higher order accurate differencing stencil does not have 
any points on the downwind side, it requires three or 
more layers of guardcells to accommodate points on the 
upwind side, which is expensive memory-wise. Thus, al- 
though higher order spatial differencing should reduce 
reflection from refinement boundaries, it is not always 
practical. 

Use of a 2nd order accurate MAD stencil, as in 
Eq. (28), for advection, avoids the above difficulties. As 
we generally locate the punctures within the finest grid 
regions, and the MAD stencil automatically reverts to 
conventional 2nd order upwinding in this region, the ad- 
vection derivative does not take any points on the down- 
wind side in the vicinity of the puncture and we find that 
stability is maintained. 

In this way, we have implemented a 2nd order accu- 
rate MAD stencil for advection. For all non-advection 
derivatives, we have implemented a linear combination of 
2nd order centered and 4th order centered stencils, as in 
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FIG. 5: Comparison of the Hamiltoniain constraint H in the 
case of 2nd order accurate, conventional differencing stencils 
and 2nd order accurate, MAD stencils. There is a single punc- 
ture black hole centered at the origin, and refinement bound- 
aries at (sc* | — 1 M y 2M, and 4 M. 

Eq. (27). As a result, reflection from refinement bound- 
aries has been dramatically reduced in the case of a single 
puncture (Fig. 5). We have found similar improvement 
in the binary black hole case as well. 

VII. CONCLUSIONS 

We have studied a problem peculiar to advection across 
a mesh interface, in the form of reflections that propagate 


at the speed {3 of the advection. Effective blue-shifting of 
the error makes it particularly egregious. We have suc- 
cessfully modeled the problem analytically, and thereby 
found a solution. Our proposed solution, that of mak- 
ing the differencing error continuous across mesh inter- 
faces, may have wide application. In addition to grids of 
non-uniform refinement, mesh-adapted differencing may 
also be appropriate for grids with multiple coordinate 
patches, where discontinous differencing error may also 
be expected. 
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